Embedding theory for excited states with inclusion of self-consistent 

environment screening 



Johannes Lischner and T.A. Arias 
Laboratory of Atomic and Solid State Physics, 
Cornell University, Ithaca, New York 14-853 

Dominika Zgid and Garnet Kin-Lie Chan 

Department of Chemistry and Chemical Biology, 
Cornell University, Ithaca, New York 14-853 

We present a general embedding theory of electronic excitations of a relatively small, local- 
ized system in contact with an extended, chemically complex environment. We demonstrate 
how to include the screening response of the environment into highly accurate electronic 
structure calculation of the localized system by means of an effective interaction between 
the electrons, which contains only screening processes occurring in the environment. For 
the common case of a localized system which constitutes an inhomogeneity in an otherwise 
homogeneous system, such as a defect in a crystal, we show how matrix elements of the 
environment-screened interaction can be calculated from density-functional calculations of 
the homogeneous system only. We apply our embedding theory to the calculation of excita- 
tion energies in crystalline ethylene. 

I. INTRODUCTION 

Neutral electronic excited states play an important role in the study of condensed matter sys- 
tems. They are probed in spectroscopic measurements, such as absorption, reflectivity or photolu- 
minescence. In addition, they are important in many technical applications, such as photovoltaics, 
laser technology or light-emitting diodes. It is therefore necessary to develop a theoretical under- 
standing of neutral electronic excitations. 

Various theoretical approaches to calculating excited state properties have been developed in 
the past, each of which usually is applicable to certain classes of physical systems. Highly accurate 
quantum chemistry methods, such as configuration interaction or coupled cluster theory, can only 
be applied to small systems containing few electrons, such as atoms or small molecules. For larger 
molecules and clusters, time-dependent density-functional theory [1, 2] yields a reliable description 
of excited states. Many-body perturbation theory, where neutral excitation energies and oscillator 
strengths are obtained by solving the Bethe-Salpeter equation for the two-particle Green function, 
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constitutes a standard tool for neutral excited states in extended periodic systems, but has also 
been applied successfully to localized systems, such as molecules [3-5]. However, many-body 
perturbation theory, while being applicable to much larger systems than the quantum chemistry 
methods, still scales quite unfavorably with system size and often further approximation, such as 
model dielectric functions, are introduced when studying chemically complex crystals [6, 7]. 

Often, in extended systems one can identify a subsystem, which is of particular interest: for 
example, an adsorbed molecule near a solid surface or a defect in a crystal. In such systems, 
the electronic excitation is often localized on the special subsystem and a high level of accuracy is 
desired for the description of the excitation. However, while standard approaches for large systems, 
such as time-dependent density-functional theory, often do not yield the required accuracy, highly 
accurate quantum chemistry approaches can only be employed when the environment, which screens 
the potential created by the localized excitation, is ignored. Embedding theories, which attempt 
to combine a high-accuracy treatment of the subsystem with a more approximate treatment of 
the environment, attempt to overcome this difficulty. Those theories have a long history, but have 
focused mostly on ground-state properties [8]. 

Regarding excited state embedding, Whitten and coworkers [9] construct an embedding scheme, 
where a configuration interaction calculation of the subsystem is based on a Hartree-Fock calcula- 
tion of a surrounding cluster. However, the description of the environment by a cluster is a severe 
approximation, which neglects long range screening effects. Carter and coworkers [10-14] building 
on previous work[15, 16] embed a configuration interaction calculation into an environment de- 
scribed by density-functional theory. Starting from a formally exact set of equations, Carter and 
coworkers [10] compute excitation energies of the subsystem by introducing an additional "exter- 
nal" potential due to the environment into their quantum chemistry calculation. However, they 
assume that the electron density of the environment in the excited state remains the same as in 
the ground state, neglecting the rearrangement of the environment electrons due to the excitation 
of the subsystem. 

In this work, we present a theory of electronic excitations of a localized system in contact with 
an extended, chemically complex environment taking proper account of the environment response 
to the excitation on the localized system. In particular, we include the environment response 
through a special screened interaction, which contains only screening processes of the environment, 
acting between electrons of the subsystem. This environment-screened interaction is obtained by 
first carrying out density-functional calculations of the homogeneous system, e.g. the defect-free 
crystal, and then subtracting out screening processes due to the region which is replaced by the 
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explicit system (e.g. defect). 

Once obtained, the environment-screened interaction can be used in various electronic structure 
methods for the localized, "explicit" system: we demonstrate how to include the environment 
response into a Green function calculation of the subsystem or into a wavefunction calculation. 

This paper is organized as follows: in Section 1, we introduce the concept of an environment- 
screened interaction and describe its use in an embedding procedure, where the explicit region is 
treated by Green function methods. In Section 2, we describe the application of our embedding 
approach to wavefunctions methods, which are typically employed in quantum chemistry calcula- 
tions. In particular, we introduce a self-consistent equation for the environment-screened potential 
caused by an external charge distribution and consequently transform this equation into a matrix 
form. We also discuss numerical difficulties of our embedding approach which are caused by (i) 
the divergence of the screened interaction and (ii) the extreme narrowness of the basis functions 
used in quantum chemistry. The first problem is solved by introducing a new integration scheme 
for integrals over the Brillouin zone, which avoids the approximations made in similar approaches. 
The second difficulty is overcome by computing the changes of the interaction matrix elements due 
to the presence of the environment, which can be evaluated efficiently in Fourier space because 
only the smooth parts of the basis functions are screened. In Section 3, we apply our theory to 
crystalline ethylene and discuss the effects of the crystalline environment on the excited states of 
these materials. Section 4 offers discussions and outlook. 

II. EMBEDDING THEORY AND APPLICATION GREEN FUNCTION METHODS 

The identification of a localized subsystem, on which the modeling effort is concentrated, marks 
the starting point in every embedding approach. In our embedding approach for excited states, 
the physical system under consideration is divided into two regions: a localized "explicit" region, 
where the excitation occurs, and the rest of the system, the environment, which screens the potential 
created by the excitation. Having identified the relevant subsystem, we employ highly accurate, 
but computationally demanding electronic structure methods, such as the Green function methods 
or wavefunction methods, to describe the localized excitation. To include the influence of the 
environment into these calculations, we modify the interaction between the electrons to contain 
all screening processes occurring in the environment. This environment-screened interaction is 
obtained from less demanding many-body Green's function methods, evaluated eventually with 
the Kohn-Sham orbitals from density-functional theory calculations. 
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A. Environment-screened interaction 

The environment-screened interaction between electrons of the localized system, where the 
excitation occurs, plays a crucial role in our embedding approach. Quite often, it is possible 
and advantageous to replace a complicated calculation of a self-consistent environment response 
by the use of an effective or screened interaction, the simplest case being the interaction of two 
point charges in a linear dielectric medium. The notion of a screened interaction also plays an 
important role in quantum many-body theory, where Feynman diagrams are used to organize 
and gain intuition about the plethora of possible screening processes. Here we also employ a 
diagrammatic approach to derive an explicit expressions for the environment-screened interaction. 

The rules for converting Feynman diagrams into algebraic expressions can be found in many 
textbooks, e.g. Ref. [17, 18]. In the diagrams below, dotted lines denote bare Coulomb interactions, 
while the fully screened (containing screening processes from the environment and the localized 
subsystem) interaction is represented by double-dotted lines. As usual, electron (and holes) are 
denoted by regular lines. The screened interaction is given by 

- ----- • O • O O • CD + 

By introducing the irreducible polarizability, represented below by a dashed square, we can 
express the infinite sum by a self-consistent Dyson equation. The irreducible polarizability contains 
all diagrams which cannot be separated into two parts by cutting a single dotted line. Dyson's 
equation for the screened interaction is given by 



= + * • (2) 

Of course, the exact form of the irreducible polarizability is unknown and must be approxi- 
mated. A well-known approximation, which has proven extremely useful in many applications of 
electronic structure theory, is the random-phase approximation (RPA) , where only simple electron- 
hole bubbles are inserted into the dotted lines in Eq. (1). 

The electron and hole of each bubble are either part of the explicit subsystem or part of the 
environment. In the following, Green function lines corresponding to subsystem electrons (and 
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holes) are colored red, environment electrons (and holes) are colored blue. Consequently, a diagram 
representing an electron-hole bubble gives rise to four colored diagrams 

____Q___ _ 

(3) 

In the first colored Feynman diagram an electron- hole pair is created in the environment, while 
in the second diagram the pair is created in the explicit region. In the third and fourth diagrams 
one particle is created in the environment, while the other is created in the explicit system. The 
contribution to the total screening from such diagrams, where environment particle lines and 
system particle lines connect, is small if there is little wavefunction overlap between subsystem and 
environment. Algebraically, such an electron-hole bubble corresponds to the product of a Green 
function G sys (r,r') describing particles in the explicit system and a Green function G env (r,r') 
describing environment particles. Because the bare Green function is proportional ip{r)ip*{r'), the 
product will vanish if environment and system wavefunctions do not overlap. 

Within the RPA, the fully screened interaction can therefore be expressed as 

= ----- +— o- +— a- +--0-0 +-0-0- +- 

(4) 

We now define the environment-screened interaction, which is represented by a blue double- 
dotted line and contains only environment bubbles, according to 

..... = +— O- +-00- +-0~0~0~ +...(«) 

Note that the fully (within the RPA) screened interaction can be obtained from the environment- 
screened one via 

..... ====== • O + --0-0 + --00-0- +■■■ 

(6) 
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This simple result has important consequences: it allows us to rigorously replace the whole 
system, consisting of the localized subsystem and the environment, by an "effective" subsystem, 
in which electrons interact through the environment-screened interaction. 

We note that this result is not limited to the random-phase approximation of the screened 
interaction, but is found quite generally whenever the irreducible polarizability is a sum of two 
pieces: one containing only environment electrons, the other only system electrons. This opens 
up the possibility of concentrating the modelling effort on the localized subsystem by using more 
sophisticated approaches for this region than the RPA. In our approach, we employ the RPA for 
the environment-screened interactions and use highly accurate electronic structure methods, such 
as Green function theory or explicit wavefunction approaches, for the localized subsystem. 

B. Embedding in Green function methods 

In this section, we demonstrate how to use the environment-screened interaction to construct 
an embedding theory, where the excitation on the localized system is described by Green function 
methods. 

The Bethe-Salpeter equation of many-body perturbation theory solves for the dressed particle- 
hole Green function [2, 3], whose poles give the excitation energies of the system, and has been 
applied successfully to extended periodic systems, but also clusters and molecules [3, 4]. 

Diagrammatically, the Bethe-Salpeter equation is given by 

where the term in parenthesis represents the Bethe-Salpeter approximation to the self-energy 
of the particle-hole Green function and contains an unscreened exchange interaction and a fully 
screened direct interaction. Typically, the fully screened interaction is obtained from the static 
dielectric matrix within the random-phase approximation [2, 3]. 

To obtain an embedding theory within the Bethe-Salpeter framework, we again color the electron 
lines red for system electrons and blue for environment electrons. Assuming that the exciton is 
localized in the explicit region, we impose that the incoming and outgoing lines in Eq. (7) are red 
system lines. The direct interaction already is fully screened containing the screening contributions 
from both the localized system and the environment. However, expressing the self-consistent Bethe- 
Salpeter equation (Eq. (7)) by an infinite sum of diagrams, we observe that the exchange interaction 
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in the two-particle self-energy gives rise to diagrams like 




(8) 



To capture such an exchange-type process, where a particle-hole pair of the localized system anni- 
hilates and re-emerges after creating an intermediate electron- hole pair in the environment, in our 
embedding approach, we replace the bare exchange diagram in Eq. (7) by an environment-screened 
exchange. 

It is easy to verify that the following Bethe-Salpeter equation for the embedded system 




indeed gives rise to diagrams like Eq. (8). 

However, apart from the aforementioned embedding approximations (neglect of diagrams with 
environment lines connecting into system lines and the assumption that the exciton resides on the 
explicit system), the proposed embedded Bethe-Salpeter equation contains one further approxima- 
tion: all environment bubbles are empty, i.e. diagrams like 



which can be constructed from Eq. (7), are absent in Eq. (9). In these diagrams, additional 
interactions between electron-hole pairs of the environment are included. Nevertheless, we expect 
that Eq. (9) captures the vast majority of all screening processes. 

Computationally, the embedded Bethe-Salpeter equation has the advantage that, once the 
environment-screened interaction is computed, Eq. (9) can be conveniently solved in a localized 
basis, such as gaussians [3], with relatively few basis functions. 



Having introduced the general notion of an environment-screened interaction and having seen 
its use in an embedded Bethe-Salpeter equation, we now move on to embed a localized system 
described by wavefunction methods, which are typically used in quantum chemistry calculations, 
into an environment described by the random phase approximation of many body theory. 




(10) 



III. EMBEDDING IN WAVEFUNCTION METHODS 
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We use the Lehmann representation of the particle-hole Green function [19] to connect the 
many-body wavefunctions and the corresponding eigenenergies to the diagrammatic Green function 
formulation of the last section. Thus, we may view the result of a full configuration interaction 
calculation as the solution of a generalized Bethe-Salpeter equation with the exact two-particle 
self-energy for the electrons of the localized subsystem. Diagrammatically, this corresponds to 




(11) 

where the term in brackets now contains more complicated irreducible diagrams. Also, the fully 
screened interaction is not given by Eq. (6) any more, but contains more complicated screening 
processes due to the exact treatment of the subsystem 

- + --0- + -0"0" + --G>- + - 

(12) 

Again, the use of the environment-screened interaction in wavefunction methods allows for the 
application of these methods to large systems containing many electrons, which traditionally are 
out of reach for these methods. In addition to the environment-screened interaction, which only 
describes response properties of the environment, we note that it is also necessary to include an 
additional external potential into the quantum chemistry calculation, which is caused by the ground 
state configuration of the environment. 

A. Self- Consistent Equations for Environment-Screened Potential 

Having demonstrated the usefulness of an environment-screened interaction for excited state 
embedding approaches, we now turn to the task of computing it. While the calculation of the fully 
screened interaction is relatively straightforward in localized or periodic systems, the calculation 
of the interaction with only environment screening proves more difficult if the localized system 
constitutes an inhomogeneity in an otherwise homogeneous system, such as a defect in a crystal. 
For this case, we demonstrate how matrix elements of the environment-screened interaction can 
be calculated from random phase approximation calculations of only the homogeneous system 
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(crystal without defect) by subtracting out screening contributions from the region where the 
inhomogeneity (defect) resides. The desired matrix elements are computed from the environment- 
screened potential resulting from the "external charge density" of a product of basis functions. This 
potential obeys a self-consistent equation introduced in this section, which - in the next section - 
is then transformed into a matrix form. 

After choosing an explicit set of basis functions gi (assumed to be real) for the quantum chem- 
istry calculation of the embedded localized subsystem, the interaction v(r,r') between electrons 
appears in the resulting matrix form of Schroedinger's equation only in interaction matrix elements 
of the form 

(Si,S2Mff3,S4> = J d 3 r J d 3 r' gi (r)g2(r)v(r,r')g 3 (r')g^r'). (13) 

In our embedding approach, we replace the bare Coulomb potential $>(r) = 

/ d 3 r'<73(r')<74(r')/|r — r'| due to the "charge distribution" 53(74 in Eq. (13) by the environment- 
screened potential <E> according to 

(5i,52|$) -> (91, 92^) ■ (14) 

$ takes into account the response of all environment molecules to the "external charge density" 
Pext = 9394 and solves the self-consistent equation 

4> = K(p ext + Xenv$), (15) 

where K denotes the bare Coulomb operator [_Kp](r) = f d s r'p(r')/\r — r'| and Xenv = Xcrys ~Xsys 
is the environment polarizability given by the difference of the polarizability of the full crystal 
(without the inhomogeneity) and the polarizability of the local subsystem. We note that Xenv® is 
the self-consistent induced charge density in the environment caused by the external charge p ex t- 
Solution of Eq. (15) is complicated by the presence of both quantities that describe extended 
systems and are best represented in Fourier space, such Xcrys, and quantities which describe the 
explicit local system and are best represented by localized functions, such as Xsys- We overcome 
this difficulty by replacing solution of Eq. (15) by a two step process: In the first step, we compute 
the effective field due to an external charge ptot allowing the whole crystal (and not only the 
environment region) to screen the external charge. Consequently, $ obeys 

& = K(ptot+Xcrvs&)- (16) 

Expressing all quantities in plane waves, this equation can be solved straightforwardly yielding 

$ = (K^ 1 - Xcrys)" 1 Ptot = K cry s p tot , (17) 
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where K crys denotes the fully screened interaction of the homogeneous system (e.g. perfect crystal 
without defect). 

In the second step, we impose that the total external charge in Eq. (16) must be given by 

Plot = Pext + Ap = p ext - Xsys®, (18) 

which involves only the localized system response and can be solved in a localized representation. 
Inserting Eq. (18) into Eq. (16) clearly gives back Eq. (15). Alternatively, we can substitute Eq. (17) 
into Eq. (18) yielding 

Ap = XsysKcrysPtot = ~Xsys^crys{Pext ~l~ Ap), (19) 

which now constitutes a self-consistent equation for Ap. In the following, we will make use of this 
latter equation. 



B. Transformation to matrix equation 

We now show how to transform Eq. (19) for the charge density Ap induced in the environment 
into a matrix equation. 

We employ the random-phase approximation for the polarizability of the explicit system 

jk 

where /; = or 1 are occupation factors, 4>i(r) and ej denote single-particle wavefunctions and 
orbital energies of the local system, respectively. From now on, we will work in the static approxi- 
mation setting oj = 0, which is well-justified for many applications [20]. 

Introducing the standard transition-space notation, where the transition of an electron from 
orbital k to orbital j is labeled by p = (k,j), we express Xsys as 

XwM= E P^)x»p;(r'), (21) 

H={v,c) 

with Xn = 4/( e j — pfx{r) = </>jJl(r)0j(r) and the indices v and c run over occupied and empty 
orbitals, respectively. Inserting this expression for Xsys into Eq. (19) yields 

Ap(r) = -J>^(r)^| d 3 r'j d 3 r" P ;(r')K crys (r' , r")\p ext (r") + Ap(r")]. (22) 
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The ansatz Ap(r) = A^p^r) then leads to the following matrix equation for the coefficients 

\ = -Xn (a» + M^A}j , (23) 

where we defined 

a„ = J d 3 r J d 3 r'p^(r)K crys (r,r')p ext (r'), (24) 
M„ v = Jd 3 rJ d\'pl(v)K crys {v,v')p u (v'). (25) 

Solving Eq. (23) we arrive at 

A = -[ X - 1 + M}- 1 a, (26) 

with xj^l = X^^iiv ■ In sum, we solved Eq. (19) for Ap in terms of the transition-space matrix M 
and the vector a, which involve the fully screened interaction K crys of the corresponding homo- 
geneous system. In the next section, we describe how these quantities are computed within the 
random-phase approximation. 

C. Calculation of M and a 

Exploiting the translational symmetry of the crystal without the inhomogeneity (i.e. the local- 
ized subsystem), we express K crys , p ext and p^ as integrals over the Brillouin zone. For example, 
is given by 

P,(r) = [ pLp^)e^\ (27) 
Jbz vbz 

where Vbz = (2vr) 3 /V ce u denotes the volume of the Brillouin zone with V ce u being the unit cell 
volume. We note that p M , q (r) = J2g /V,q(G) exp(iG • r), where G denotes a reciprocal lattice 
vector, is a lattice-periodic function. 

To compute K crys [Eq. (17)], the polarizability of the crystal is taken in the static random-phase 
approximation [21, 22] and explicitly given by 

tc^ n>\ 4 \^ kj e ^(q+ G )- r jc, k + q) (c, k + q|e ^ (q+ G, )^ r, \v, k) 

Xcrys,^, I* J - ^ ~ " ' 

q cell tX ^,k-e c , k +q 
where N q is the number of sample points in the Brillouin zone and the indices c and v run over 
conduction and valence states, respectively. In the molecular crystals we are interested in orbitals 
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on neighboring molecules overlap very little, allowing for additional simplifications: (i) we assume 
flat bands, i.e. e„ 5 k = e n , and (ii) we assume that the crystalline Bloch states can be constructed 
from the wavefunction at the T-point through a tight-binding procedure via 

nk (r) = ^=Y1 ^ V*(r " R )> ( 29 ) 

where R runs over all N q units cells comprising the crystal and w n (r) = ^nk=o( r ) m the unit 
cell surrounding the origin and w n (r) = in all other cells. Formally, Eq. (29) describes 
the Wannier transformation from localized orbitals or Wannier functions to Bloch states us- 
ing an approximate Wannier function w n [note that the exact Wannier function is given by 
w n R, = \/y/W q x J] k exp(-ik • R)V>nk]- 

The tight-binding ansatz eliminates the summation over the Brillouin zone in Eq. (28) and 
yields 

X<WG,G') = ^^He-^+ G )- r |c)^— (c|e^ +G ')- r >), (30) 
v ceii ^ e v -e c 

where {v\e~ t( - Ci+G ^ r \c) = J v d 3 r0*(r)e~ l ( q+G )' r ^ c (r). We point out that the tight-binding ap- 
proximation for the crystal polarizability is invoked only for numerical convenience and that all 
following expressions can also be evaluated using the full RPA polarizability, given in Eq. (28), 
instead. 

Finally, we express K crys in terms of the symmetrized dielectric matrix e q (G, G') = 5gc — 
< /2 (G) Xcn/s , q (G,G')< /2 (G') via 

K cryS)(i (G, G') = ^ q /2 (G)e- 1 (G, G')K^\G% (31) 

with Kq{G) = 47r/|q + G| 2 . The components of M and a can now be expressed completely in 
Fourier space and are given by 

a M = V cM / ^ E />; q (G)iW, q (G, G')fe t , q (G'), (32) 
Mfj, v = Veen [P-Y, ^, q (G)iv cr , s , q (G, G')^, q (G'). (33) 

J BZ GG' 
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D. Environment-screened interaction matrix element 



The final expression for the environment-screened matrix element is 



(9i, 92 |$) 



= {gi,92\K crys \p ext (3,4:) + Ap(3,4)} 

= (9i,92\K cr y S \g 3 ,g 4 } + {gi,g 2 \K crys \^A IJi {?>,4)p IJ ) 



= (9i,92\K crys \gs,g4) - a'(l, 2)Xa(3, 4), 



(34) 



where we defined X = [x" 1 + M]~ l and we used Eqs. (17), (18) and (26). We also now explicitly 
write a(3, 4) to denote explicitly the dependence on the basis functions gi (similarly for Ap and 



In this section, we discuss the two major numerical difficulties of the described embedding 
approach: (i) the integrations over the Brillouin zone in Eqs. (32) and (33) and (ii) the handling 
of extremely narrow basis functions. 

Regarding the first issue, we note that the Brillouin zone integrals in Eqs. (32) and (33) must 
be handled with great care, because certain elements of K crys ^ diverge as |q| approaches zero: 
the head K crySt(l (0, 0) diverges as 1/q 2 and the wings, K crySj q(0,G') and K cryStCl (G, 0), diverge as 
l/|q|. Although these singularities are integrable, they can cause extremely slow convergence when 
the integral is approximated by a discrete sum over sample points in the Brillouin zone. 

Most approaches [23, 24] use a regular mesh including the origin to carry out the Brillouin zone 
integration. The contribution from the volume element at the origin is obtained by approximating 
the smooth part of the integrand by its value at the origin, replacing the volume element by a 
sphere of equal volume and finally analytically integrating the divergent part of the integrand in 
the spherical volume. While numerically efficient, we expect the error associated with this scheme 
to converge quite slowly. Also, in anisotropic systems the smooth part of the integrand might not 
even be well-defined at the origin and instead vary strongly depending on the direction of approach 
towards the origin. 

In contrast, in our approach we first isolate the singular parts of the integrand by splitting the 
integral over the Brillouin zone into two parts. The part containing the singularity is transformed 
to spherical coordinates thus removing the singularity, while the other part is evaluated using a 



A). 



E. Numerical challenges 
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regular grid over the Brillouin zone. For a generic function /(q), we write 




(35) 



where a is a parameter which can be adjusted to optimize convergence. Note that the integrand 
in the first term on the right hand side is well-behaved at the origin, even if /(q) diverges as 1/q 2 . 
We use a simple regular grid over the whole Brillouin zone to evaluate this term. The integrand 
of the second term is divergent at the origin if /(q) diverges, but vanishes rapidly as |q| increases 
because of the fast decay of exp(— a 4 q 4 ). The singularity can be removed by transforming the 
integral to spherical coordinates yielding 



where the integrand in parentheses is now well-behaved at the origin and relatively smooth, such 
that we can use a Gaussian quadrature scheme for accurate evaluation. We also use Gaussian 
quadrature to carry out the integral over spherical angles. 

Our scheme thus avoids the approximations described above. In particular, the contribution 
from the singular region surrounding the origin is computed with high accuracy. Also, the smooth 
part of the integrand is never evaluated at the origin, where it is not necessarily uniquely defined, 
making our scheme ideal for the study of anisotropic systems. 

Regarding the second issue raised above, the difficulty of handling extremely narrow basis 
functions, we point out that, in actual density-functional calculations, the Kohn-Sham orbitals are 
represented on discrete real-space or Fourier space grids, which are typically too coarse to faithfully 
represent the narrow basis functions used in quantum chemistry calculations. However, because 
the crystal polarizability [Eq. (28)] entering the fully screened interaction is computed from Kohn- 
Sham orbitals only Fourier components of the "external charge density" (caused by a pair of basis 
functions) which belong to the grid of the density-functional calculation are screened, suggesting 
that the crystal-screened interaction can be decomposed as 



where AK crys now only acts on the space of functions representable on the grid of the density- 
functional calculation. Instead of the full environment-screened interaction matrix element 
[Eq. (34)], the calculation of which would require grids fine enough to represent the narrow basis 
functions, we compute the difference (51(72 1 AK crys \ g% g^) , which allows us to work with the grids 
of the density-functional calculation only. This difference has to be added to the bare Coulomb 
matrix element, which is computed with high accuracy in existing quantum chemistry codes. 




(36) 



K, 



crys 



= K + AK ( 



crys j 



(37) 



15 



F. External potential due to environment 

In addition to the environment-screened interaction, we also include a static external potential 
caused by the ground state configuration of the environment into the quantum chemistry calcu- 
lation. The inclusion of this additional external potential is necessary, because the environment- 
screened interaction only accounts for response effects in the environment, but not for the static 
effect of the environment ground state configuration on the electrons of the localized subsystem. 
The proper definition of the external potential requires great care in order to avoid "double count- 
ing" of interaction processes between the explicit system and the environment. 

To reproduce the correct change in electron density Sp in the localized system system upon 
embedding, the external potential (p ex t must fulfill 

5p = Xscr4>ext, (38) 

where Xscr denotes the response function of the localized subsystem with an environment-screened 
interaction between the electrons and we assumed that 4> ex t is relatively weak. If the wavefunction 
overlap between subsystem and environment is small, we can relate <f> ex t to the total electrostatic 
potential 4> env due to the environment ground state configuration. cj) env is given by 

(frenv = Kpenvi (39) 

where the environment charge density p env is the sum of electronic and nuclear contributions. From 
its definition, (f) env produces the correct change in density Sp, when applied to the subsystem with 
unscreened interactions between the electrons. If the potential is weak, this implies 

bp = XO 4>env, (40) 

where xo denotes the response function of the subsystem with unscreened interactions. Combining 
Eq. (38) and Eq. (40), we find 

4>ext = XscrXo4>env, (41) 

which allows for the calculation of the correct external potential from the electrostatic potential 
caused by the environment ground state configuration. 

IV. APPLICATION TO ORGANIC CRYSTALS 

In this section, we apply the described embedding approach to study electronic excitations in 
crystals consisting of short organic molecules. In these crystals, we take the molecules of a single 
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unit cell as the localized subsystem and all other molecules as the environment. 

As a test case, we investigate crystalline ethylene. Crystalline ethylene being one of the simplest 
organic molecular crystals has been studied for a long time[25-27]. 

The excitations of isolated short organic molecules have attracted considerable attention because 
of the non-trivial character of their low-lying excited states [28-30]. For example, contrary to naive 
expectations the lowest excited state of the octatetraene molecule is not a one-electron (HOMO- 
LUMO) excitation, but has a HOM0 2 -LUM0 2 double-excitation character [29, 30]. 

Upon embedding the molecule in a crystalline environment, we expect the molecular excitations 
to be altered because of (i) the screening by the crystalline environment and (ii) the derealization of 
the excited state over several neighboring molecules. For crystals of short acenes, Ambrosch-Draxl 
and coworkers [31] found that there is very little derealization of the excitation. We therefore only 
include the two inequivalent molecules of a single unit cell into our embedded quantum-chemistry 
calculation. 

A. Computational details 

We carry out density-functional calculations of ethylene crystals using the generalized gradient 
approximation [32]. We employ Kleinman-Bylander pseudopotentials [33], a plane wave cutoff of 
40 Hartree and a 2 x 2 x 2 Monckhorst-Pack [34] kpoint grid. For ethylene, the crystal structure is 
well-known from x-ray and neutron diffraction and we employ the lattice parameters and atomic 
positions provided by Ref. [35] in our calculation. We then relax the ionic positions and compute a 
large number of empty bands, which are needed for the calculation of the polarizability [Eq. (20)]. 

Next, we evaluate the environment-screened interaction matrix elements, Eq. (33). In this calcu- 
lation, we only employ a subset of the reciprocal lattice used in the density-functional calculations 
consisting of ~ 1000 reciprocal lattice vectors. We find that the matrix elements are sufficiently 
converged if 200 empty bands are used. For the Brillouin zone summation, we employ a regular 
grid containing 512 kpoints for the first integral in Eq. (36) and 270 kpoints for the second integral, 
which contains the singularity. Finally, we approximate the external potential by the electrostatic 
potential of the environment, (p ex t ~ 4> e nv 

To compute excited states of the localized subsystem, we carry out configuration interaction 
calculations using the 6-3 1G basis set. 
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bare 


screened 




[eV] 


[eV] 


state 1 


17.432 


17.384 


state 2 


17.436 


17.554 


state 3 


19.004 


19.004 



TABLE I: Lowest excitations energies (measured from the ground state energy) in crystalline ethylene. The 
first column contains excitation energies obtained from configuration interaction calculations of the two 
molecules in a unit cell without any environment effects. The second column contains excitation energies 
for the unit cell embedded in the crystalline environment. 

B. Results 

Table I shows our results for the energies of lowest lying excitations in crystalline ethylene. 
We observe that the changes due to the environment screening are smaller than one percent. 
Interestingly, we find that the inclusion of environment screening lowers the energy of the lowest 
excited state, but increases the energy of the second lowest excited state. The energy of the third 
lowest excited state remains constant. 

C. Discussion and Outlook 

In this work, we have described a general framework for embedding theories of excited states 
in physical systems where the excitation is localized on a small subsystem. Contrary to previous 
approaches, we fully include the self-consistent screening response of the environment by using an 
effective interaction between electrons in the explicit subsystem. Once obtained, the environment- 
screened interaction may be employed in various highly accurate electronic structure description of 
the localized subsystem: we demonstrate how environment effects can be incorporated into Green 
function methods or wavefunction methods, which are typically employed in quantum chemistry. 

We apply our embedding theory to the calculation of excited state energies of crystalline ethylene 
and find encouraging results. 
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